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ABSTRACT 

A simple, yet general, formalism for the optimized linear combination of astrophysical images is constructed 
and demonstrated. The formalism allows the user to combine multiple undersampled images to provide over- 
sampled output at high precision. The proposed method is general and may be used for any configuration of 
input pixels and point spread function; it also provides the noise covariance in the output image along with 
a powerful metric for describing undesired distortion to the image convolution kernel. The method explicitly 
provides knowledge and control of the inevitable compromise between noise and fidelity in the output image. 

We present a first prototype implementation of the method, outlining steps taken to generate an efficient 
algorithm. This implementation is then put to practical use in reconstructing fully-sampled output images 
using simulated, undersampled input exposures that are designed to mimic the proposed Wide Field InfraRed 
Survey Telescope (WFIRST). We examine results using randomly rotated and dithered input images, while 
also assessing better-known "ideal" dither patterns: comparing results we illustrate the use of the method as a 
survey design tool. Finally, we use the method to test the robustness of linear image combination when subject 
to practical realities such as missing input pixels and focal plane plate scale variations. 

Subject headings: dark energy — methods: data analysis — techniques: image processing — techniques: 
photometric — weak gravitational lensing 



1. INTRODUCTION 



Weak gravitational lensing (e.g. I Schneide^l2006l l) is a pow- 
erful means by which to probe the distribution of matter in the 
Universe. As such, it has been identified as among the most 
promising methods by which to constrain the growth rate of 
matter structure in the Universe, and thereby to test models 
of evolving dark energy and modifications to general relativ- 
ity (see, e.g.. lAlbrecht etai]|2006l l2009t iPeacock et al.ll2006i 
and references therein). However, weak lensing is arguably 
the most technically challenging of the cosmological probes 
from an image analysis standpoint, due to its extreme sensi- 
tivity to any systematic errors made when recovering galaxy 
shape information. Much work has gone into understanding 
how to accurately recover tiny, coherent gravitational shear 
signals from large numbers of noisy galaxy images that have 
been convolved with the instrumental Point Spread Function 
(PSF) (Hey mans et a l. 2006; Massev et al. 2007; Bridle et al] 
[2010; Kitc hing etalJi2010 ). 

In the analysis of astronomical images, ensuring the ad- 
equate spatial sampling of data by pixels of finite size and 
spacing is a key concern. Ideally, images should be sampled 
at or above the Nyquist-Shannon sampling rate for th e band 
limit set by the optical response of the system (see e.g. Marks 
120091) . so that the full continuous image can be determined 
from the discrete pixel samples. Throughout this Paper we 
will define a single plane wave to be cx e^'^'" so that u has 
units of cycles per arc second. In general, if an image contains 
only Fourier modes whose spatial frequency u is no larger 
in magnitude than Umax then the Nyquist criterion demands 
sample spacing P satisfying P < l/(2Mniax)- An image sam- 

|browe@caltech.edu; chirata@tapir.caltech.edu; jason.d.rhodes@jpl.nasa.gov 
' Jet Propulsion Laboratory, California Institute of Technology, 4800 

Oak Grove Drive, Pasadena, CA, 9 11 09 

- California Institute of Technology, 1200 E California Boulevard, 

Pasadena, C A, 91106 

Department of Astrophysics, Califoniia Institute of Technology, M/C 

350-17, 1216 E California Boulevard, Pasadena, CA, 91106 



pled more finely than the critical rate 1 / (2M,„ax) is referred to 
as oversampled (and one sampled at the critical rate is criti- 
cally sampled). Since it is possible to reconstruct the entire 
original image from oversampled data, operations such as in- 
terpolation/regridding, rotation, and translation can be carried 
out with no pixelization artifacts. This makes oversampled 
data the preferred input for most precision image analysis ap- 
plications, including weak lensing. 

The ideal location for weak lensing observations is a space- 
based telescope, where one is free from the blurring effects 
of the Earth's atmosphere and can achieve a level of sta- 
bility of the optics and hence the PSF that is impossible 
from the ground. Weak lensing is thus a key project for 
proposed space-based imaging missions such as the Wide 
Field InfraRed Survey Telescope (WFIRST: Blandfor d et al.l 
120 lot the reference design is based on the Joint D ark En- 
ergy Mission/JDEM- Omega concept: lGehrelsir2010l) and Eu- 
clid (Refre gier et al.l 12010) . However, in both cases, practi- 
cal design considerations prohibit oversampling at the native 
pixel scale of the system. The optics of a space telescope de- 
liver a PSF that preserves Fourier modes out to Umax = D /\ 
where D is the outer diameter of the primary mirror and A 
is the wavelength of observation; the high spatial frequen- 
cies may be suppressed by e.g. charge diffusion, but in most 
cases Umax is still large. Indeed, preserving the high spa- 
tial frequency components of the image is a major reason to 
choose a space-based platform. If one is to oversample the 
image at the native pixel scale, then, one is forced to choose 
very small pixels. However, there are competing considera- 
tions that drive one to larger pixels, including (i) the desire 
for large field of view within engineering or cost constraints 
on the number of detectors; and (ii) the high read noise of 
near-infrared detectors, which results in increased photomet- 
ric errors as light from an object is spread over more pixels. 
For these reasons, both JDEM-Omega and Euclid chose pixels 
larger than 1 / (2umax)- Images from such systems, sampled at 
the native pixel scale, are undersampled. Undersampled im- 
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ages suffer from aliasing, in which each Fourier mode of the 
observed image contains contributions from several Fourier 
modes of the original image: in this case the original image 
cannot be unambiguously reconstructed. Operations such as 
interpolation, translation by a non-integer number of pixels, 
etc. are not mathematically possible on undersampled data, 
which represents a major difficulty for image analysis. 

Although proposed space-based dark energy missions are 
unlikely to oversample their imaging data on a given single 
image, typical design concepts allow accurate pointing capa- 
bility (e.g. Gehrels 2010). What is required is thus a means 
of reconstructing oversampled images from two or more un- 
dersampled, dithered images. The simplest technique to im- 
prove sampling is the ideal subpixel dither: for example, one 
may take 4 exposures whose relative offsets are (0, 0), (i, 0), 
(0, i), and (i, i) pixels (a "2 x 2" dither pattern). With the 
use of these multiple exposures, where the observed image de- 
pends on different linear combinations of the aliased modes, 
one may use linear algebra te chniques to s eparately "solve 
out" the aliased Fourier modes (lLauerl ll999a b) so long as the 
native pixel scale P < 1/umax- The linear algebra techniques 
allow us to analyze any subpixel dither pattern: notable ex- 
amples include the "3 x 3" pattern (9 exposures); subpixel 
dithers with diagonal grids such as the "\/2 x v^" pattern 
[where the relative offsets are (0, 0) and (i, i) pixels] or the 

"\/5 X -\/5" pattern (see Section |S!2] i; or N exposures with 
arbitrary offsets. Given any such subpixel dither pattern, one 
can determine whether the linear system for the Fourier mode s 
of the input image can be solved. The output image from this 
approach has the same PSF as the input images. 

However, it is likely that WFIRST - or any wide-angle sur- 
vey that attempts to fill in gaps between detectors - will end 
up with a survey pattern that contains large slews (of order a 
detector width) between each exposure. Under such circum- 
stances, the different input images that are being used to build 
the final, oversampled output image will be sampled on dif- 
ferent grids whose phase cannot be controlled, will have dif- 
ferent field distortions, and will have "holes" due to cosmic 
rays and detector effects. Moreover, the input PSFs will be 
different even if the telescope optics are perfectly stable, be- 
cause each observation of a galaxy will be made in a different 
part of the field. The Fourier space/linear algebra techniques 
are not well-equipped to handle such situations efficiently. 
There are stable, commo nly-used image combina tion algo- 
rithms such as Drizzle (Fruchter & Hook"2002) with ma- 
ture accompanying scripts (e.g. Koekemoer et al. 2002) that 
will produce an output image for such inputs. However, these 
algorithms do not give full control of the PSF in the output 
image, or even guarantee that this PSF is constant across the 
image of a galax}|3, which should be seen as a prerequisite for 
high-precision shape measuremen t for a space weak lensing 
mission. Recently Fruchteil (|201 li) proposed a new method 
(iDrizzle: iterative Drizzle) that tackles some of these is- 
sues using an iterative application of DRIZZLE commands. It 
will be interesting in the future to compare results using this 
algorithm to those of this Paper 

A third approach to combining images would be to adapt 
methods from Cosmic Microwave Background (CMB) ex- 

* Special cases such as the 2x2 ideal subpixel dither with identical input 
PSFs may be exceptions because of special symmetiies, e.g. translation by 
half-pixels. However, our interest is in the combination of images with no 
such symmetries. 



periments. CMB experiments typically have only a few pix- 
els in their focal plane (the Wilkinson Microwave Anisotropy 
Probe AVMAP had a total of 10 differencing assemblies: e.g. 
Jarosi ket al.l 1201 ll and references therein) and use a large 
number of scans across the sky to build up an output map. 
Mapmaking for CMB experiments is usually done by taking 
the vector d of time-ordered data from a small number of de- 
tectors, and writing it in the form: 

d = Ms + 7/, (1) 

where s is a vector of temperatures in pixels on the sky, rj is 
the noise, and M is the rid x np mapping matrix (rid being 
the number of samples of data and rip the number of pixels 
in the output map). In its simplest form M is a sparse ma- 
trix containing Mij = 1 if the ith sample is collected when 
looking at the jth pixel, and otherwise, but many general- 
izations are possible (e.g. for differential experiments, meth- 
ods that mask bright or variable sourc es so they do no t con- 
taminate the data processing, etc. see iHinshaw et al.l 12007). 
Taking the known noise covariance matrix N = {rjrj^), one 
may obtain a least-squares solution for the sky map s: this is 
(MTN-iM)-iMTN-id. 

However, this method (as written here and in most of the lit- 
erature) does not trivially incorporate a variable PSF (beam), 
requires the output s to have pixels much smaller than the 
beam, and then requires at least several samples per output 
pixel - this is a luxury that is not present in optical astronomy 
with sub-arcsec PSF sizes. This problem can be solved by al- 
lowing M to contain interpolation coefficients, but an imple- 
mentation of this approach is not yet public, and is typically 
unnecessary for the arcminute beam sizes commonly seen in 
CMB experiments. 

In this Paper, we take the first steps toward developing a 
general linear algorithm that attempts to combine the best fea- 
tures of these approaches. As in the Fourier/linear algebra 
approach (Lauer 1999a b), we d esire a well-controlle d final 
PSF. As is the case for DRIZZLE (Frucht er & Hook ' 2002). we 
want to be able to use arbitrarily placed samples if possible. 
Sometimes, reconstructing a fully-sampled image with a con- 
trolled PSF is impossible given the positions of the samples 
provided; unlike DRIZZLE, we desire a method that alerts the 
user when this happens. Finally, as in the CMB approach, we 
aim to do something that minimizes output noise using some 
variant of the least-squares method. 

In Section |2] we derive a linear formalism that can meet 
these aims, and in Section [3] describe a prototype implemen- 
tation of the method called IMCOM, which is available freely 
from the authors on request. In Section |4] we illustrate the 
process of combining multiple images using IMCOM via a de- 
tailed worked example, before turning in Section |5] to exam- 
ine a set of realistic potential observing scenarios for a space- 
based dark energy mission such as WFIRST. While it should 
be stressed that the dither scenarios and PSF patterns used in 
these tests do not necessarily reflect any firm plans for the 
WFIRST mission, they are designed to approximate such a 
mission in a broad sense. In Section|6]we discuss the comput- 
ing resources necessary to process imaging data from a dark 
energy mission using the IMCOM technique, before drawing 
conclusions in Section|7] 

2. LINEAR FORMALISM 

The problem is to combine several undersampled input im- 
ages into a single oversampled image. The input images may 
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be written as a vector of intensities li of length n, where n 
is the total number of usable pixels. For example, if the in- 
puts are four 128 x 128 pixel images with no defects, then 
n = 4(128)^. The inputs are sampled at their pixel centers 
We also suppose that the PSF (including jitter, op- 
tics, and detector response) at separation s is Gi{s). That is, 
the intensities are given by 



/(r')G,(r,-r')dV + 77, 



(2) 



where rji is the noise with {rji) = and some covariance ma- 
trix A'ij = {rjirjj). The formalism can rigoroursly treat any 
general noise covariance matrix Nij, although in most cases 
Nij is close to diagonal (with small off-diagonal correlations 
due to, e.g., corrections for inter-pixel capacitance, electronic 
1// or flicker noise, etc.). The function /(r') describes the 
physical image on the sky. 

We note here that the PSF Gi (r) is assumed as being fully 
specified, in advance, at the position centered on each input 
pixel Ti, in each contributing exposure individually. This 
knowledge must be expected to come from either a well- 
motivated optical model or the image data 7^, or both. Of 
course, images of point sources (stars) within Ii will be sub- 
ject to the same degree of undersampling as the rest of the Ii 
pixels. This issue can be countered if adopting a fully para- 
metric model and fitting for Gi(r) from stellar images (e.g. 
iMa et a l. 2008). However, the use of a non-parametric model 
that potentially contains defects due to aliasing will under- 
mine the success of any image combination. The problem 
of recovering Gi (r) from undersampled stellar images is dis- 
cussed further in Section |7] 

The output is to be a synthesized, fully sampled image 
on a grid of pixel centers Rq,, where a ~ 1, . . . , m. For 
example, if the synthesized image is to be 256x256 then m = 
256^. We wish for Ha to be as close as possible to the target 
image J a, which is defined as: 

Ja= I /(r')r(R„ - r')d\'. (3) 

Here F is the desired effective PSF of the synthesized image, a 
free choice for the user A well-motivated choice for F might 
be Gi (r) itself (if it is approximately constant), preserving the 
input PSF as much as possible. This is the approach taken in 
later Sections, but in general F may be freely chosen to ad- 
ditionally filter the output image as desired (see, e.g.. Section 

Es). 

2.1. Notation 

Fourier transform pairs are written with Fourier space in 
units of cycles per arcsecond rather than radians per arcsec- 
ond. That is. 



/(u)=J-[/(r)] (u) 



/(r) 



-27riu r j2 



(4) 



for the forward Fourier transform, and 



/(r) = /(u) (r) = / /(u)e^^'"'-d2u (5) 

for the inverse Fourier transform. 

The Einstein summation convention is not used: all sum- 
mations over indices will be either be written explicitly or, 
for brevity, as implied matrix multiplications upon matrix and 



vector objects in bold face. Latin indices are used for input 
vectors and in derived objects where reference is made to in- 
put pixel locations, and Greek indices are similarly used to 
refer to output pixel locations. 

2.2. Linear Solutions 

We have defined the fully-sampled output in relation to 
the target image Jq, of equation (|3]l, to which we intend 
to be as close an approximation as possible. We then desire 
Ha to be linearly related to the input images Ii for several 
reasons. Most importantly, it ensures that the final PSF in 
Ha is independent of image brightness, so the stars that were 
not used in the PSF determination are tracers of how well the 
image synthesis worked. Linearity simplifies the analysis of 
systematics. Furthermore, it dramatically reduces the space 
of possible algorithms and allows a systematic search of the 
space that remains. We fu rther note that some pre -existing 
algorithms (e.g. DRIZZLE: iFruchter & Hookll2002h are also 
linear and hence are included in our search space. It should 
be made clear, however, that detectors usually display some 
level of non-linear behavior (e.g. non-linear gain), changing 
the PSF for bright objects. Corrections for this non-linearity 
are typically applied in image pre-processing at the raw pixel 
level, and we assume that any such non-linearity will have 
been duly corrected in Ii before starting to combine images. 

A general linear mapping from / iJ is: 



H„ 



Taili ; 



(6) 



where Tai is an m x n matrix. We now need to construct 
an appropriate objective function for the performance of the 
matrix Tai and extremize it to find the optimal linear transfor- 
mation we desire. 

2.3. Objective Function 

We note that for a matrix Tai, an error map can be con- 
structed: 

Za — Ha J a 

= f f{r')La{Ra - r')d\' +y^Ta^m 

where we have defined the leakage function La as 

i„(R„ - r') = Y,Ta^G,{v, - r') - r(R„ - r'). (7) 

i 

The leakage is simply the difference between the desired PSF 
F and its actual reconstructed counterpart in Ha ■ We will seek 
to minimize this difference. Taking equation (|7]i, we may then 
separate out the deterministic part. 



{Za) = [ /(r')ia(r'-Ra)d2r', 

JR2 



(8) 



from the stochastic contribution to the residuals as specified 
by the noise covariance. 



Cov(Zq, Z, 



^ ^ TniTi 



(9) 



We will also seek to minimize the noise variance Yiaa for each 
output pixel value Ha ■ Noisy output images will be those that 
combine input pixels from Ii with large positive and negative 
weights in Tai- Such solutions will tend to be unstable, which 
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gives a further practical reason for wishing to minimize 
where possible. 

This now allows the construction of an objective function 
W, i.e. a function to be minimized in solving for Tai- We will 
choose a function that is a simple sum over values of a, so 
that if we wish to we may write W — J2a where Wa 
depends only on ^^(s) and E^q. In other words, we choose 
to optimize each output pixel independently: one minimizes 
W by minimizing each Wa ■ 

This choice is made both for computational simplicity, and 
also so that the optimization does not mangle good regions of 
the synthesized image to try to clean up bad regions (e.g. the 
edges of frames or the occasional intersection of multiple cos- 
mic ray tracks). We desire for Wa to be a quadratic function, 
so that linear algebra methods can be used to solve for Tai, 
and for it to be translation-invariant (i.e. we care only about 
the normalization and shape of the leakage function La, not 
about its position on the sky). The most general such function 
is 



Wa — Ua + nY^aa, 



(10) 



where k > is a Lagrange multiplier, and we have defined 
the leakage objective Ua as 



Ua = 



T{s' -s)La{s)La{s')(fs(rs', (11) 



where T(s) is a real symmetric positive-definite kernel, i.e. 
T(-s) = T(s), 3T(s) = 0, and T(s) > 0. We note that 
for the simplest choice of T(s), the Dirac delta function 6(s), 
this expression becomes 



Ua 



[La{s)fd''s. 



(12) 



The leakage objective Ua is thus, in the case T(s) = S{s), 
a very simple, direct, quadratic measure of the local leakage 
at the point Rq . Returning to the general case, equation ([Til 
may also be written conveniently in Fourier space: 



Ua - 



T{u)\Lain)\^d^u. 



(13) 



It will be seen this expression for Ua is computationally con- 
venient. 

The role of the Lagrange multiplier k in equation (fTOl l mer- 
its some discussion here. It can be seen that k may be used to 
weight the relative importance of the two terms contributing 
to the objective function Wa- These relate to two key proper- 
ties of the output image: 

• The fidelity of the ensemble mean Ha to the target im- 
age of equation (O, on which Ua depends. 

• The noise variance in the output image, S]q.q. 

To ensure fidelity we naturally wish to emphasize the mini- 
mization of Ua relative to J^aa, implying k <C 1, and vice 
versa if we wish to suppress noise. Thus the minimization of 
Wa to solve for Tai, subject to the freely-chosen Lagrange 
multiplier constraint k, gives the user control of the balance 
between image fidelity and noise at each output pixel Ha ■ 

2.4. Minimization of the Objective Function 



We want to minimize Wa over the space of Tai- To do this, 
we first take the Fourier transform of equation (|7): 



Lain)-. 



X e 



27riu-R, 



-27riu-r, 



G*(u)-e- 



r*(u) 



(14) 



The leakage objective is then 



Ua 



G*(u)-e- 



r*(u) 



X t(u)d2u. 



(15) 



Recalling that Tai is real, and using equations (|9]l & (fTOl l. this 
expands to give the combined objective function: 

Wa^Y. (^"^J" + ^^N,,)Ta^Ta,+Y, Ba^Ta^+Ca, (16) 
ij i 

where the coefficients are given by the system matrices: 

Aa^,^ [ t(u)e2-"-('--'--)G,r(u)G,(u)d2u, (17) 

Ba^^-2 f t(u)e2-"("-'-i^°)G,(u)f*(u)d2u, (18) 

JR2 



Gq 



f{u)\t{u)\'d'u. 



(19) 



Here Aaij is an n x n, positive definite symmetric real ma- 
trix, Bai is a real matrix of dimension m x n, and Ca is 
a positive real scalar. (The integrals can be seen to be real 
since the integrands in Aaij, Bai, and Ca transform to their 
complex conjugates under change of the integration variable 
u — ^ — u.) We note that direct expansion of equation ( fTsT i 
would have lead to a Hermitian matrix; however, since Tai is 
real the contributions to Ua from the imaginary, antisymmet- 
ric part vanish. 

Any function written in the form of equation (fT6t can be 
immediately be minimized to give 



Tar- 2 H 



[(A„ + kN)- 



B. 



(20) 



where (A^ + kN)~^ is the n x n matrix inverse of Aaij + 
KNij . In practice it will be necessary to restrict the synthe- 
sized image pixel a to only depend on nearby input pixels i 
so as to reduce the dimensionality of the matrix Aaij ■ If this 
is done, however, Aaij actually depends on a which is not 
true if the whole image is used. Since the whole image may 
have n ^ 10^, this reduction of dimensionality is required in 
order to make the approach computationally feasible. These 
issues are discussed further in Section |3] 

We note that equation (l20l l is a general result and can be 
used in principle for any system of non-ideal dithers, rotated 
images, and even images with defects (the "lost" pixels are 
simply removed from the system matrices). Some of these 
cases will be explored in Sections|4]&|5] 

The noise variance in each output image pixel is given by 
the diagonal terms i n the full covariance of equation (|9]l. As 
discussed in Section l231 the noise properties can be adjusted 
to using the Lagrange multiplier k. In Appendix |A] we de- 
rive in detail the limiting behavior of output images across the 
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range of k values. These results make it clear that, in setting 
K, some compromise is desired between the limiting cases of 
maximum fidelity and minimum noise. In some cases, even 
if K is taken to zero (i.e. the optimization is sensitive only to 
mean output image fidelity and not to noise) the leakage func- 
tion will not approach zero. In these cases, reconstructing 
an unbiased synthesized image is not mathematically possi- 
ble using a purely linear combination of the input pixels, and 
limK^o+ Ua represents a fundamental limitation. 

The choice of k may be made globally or from pixel to 
pixel. Setting a single value of k is certainly the simplest ap- 
proach, but it is advantageous to allow the user to specify the 
properties of their input image in terms of key deliverables 
such as the noise variance Tiaa or leakage objective Ua- One 
strategy for image combination would therefore be to regard 
Eq-q as fixed (using this to set the value of k from pixel to 
pixel) and then obtain the highest-fidelity mean image consis- 
tent with this noise value. The opposing strategy would be to 
set a tolerance on the image fidelity by requiring Ua to be less 
than some threshold value, and then finding the least noisy re- 
construction consistent with this threshold. We now consider 
these strategies and give more detail about practical issues in 
the implementation of the algorithm. 

3. IMPLEMENTATION 

In this Section we describe a prototype implementation of 
the image combination method presented in Section |2l de- 
scribing some of the numerical techniques employed and is- 
sues raised in finding solutions for T^i. We have named this 
prototype package IMCOM (IMage COMbination), and plan 
to make it available to the public either in its current version 
or as part of a more developed software package at a later 
date. The code is written in FORTRAN 95, support s multi- 
thread ing via OpenMP (see, e.g.. Chap man. Jost. & PasI 
l2007h . and uses only public software libraries, some of which 
are a vailable in specific hardw are-optimized version s: CFIT- 
SIO (Penc3^9), FFTW3 (Frigo & Johns onll2005l). BL AS 
jBlack ford et al. 2002) & LAPACK ( Anders onet al.il999h . 

3. L Imcom Program Outline and Inputs 

The Imcom code solves for T^i by varying n o n an output 
pixel-by-output pixel basis as described in Section l24l so that 
K = Ka- The user may specify a maximum tolerance value of 
either the noise variance or leakage objective: these we label 
S™™ and C/™^", respectively. The algorithm may then be used 
to solve for k„ so that < S^a or Ua < C/"'% while 
simultaneously minimizing the corresponding Ua or Sqq,, re- 
spectively. This procedure leaves the user free to specify an 
appropriate compromise between noise and fidelity for the ap- 
plication in question. 

In practice it is desirable to have an iterative solution stop 
within finite time, and this can be achieved by specifying 
an acceptable range for the solutions, AE™™ and At/™". 
The algorithm can then be instructed to finish iterating once 

aa ~ ^^aa <^ ^aa ^ ^aa Or Ua ~ LlU a <■ ^ a S 

C/™^". This then allows total control over either the noise or 
fidelity of the output. 

Apart from these tolerances on the properties of the output 
image, there are other choices that the user needs to make 
before the problem is fully specified. These we refer to as 
"soft inputs", since the user has some freedom in the choice 
of their functional forms or values: 

i) The synthetic convolution kernel r(r) of the desired 



output image, 
ii) The desired output sampling locations Ha- 
ni) A choice of objective kernel T(r) for evaluating Ua- 
The simplest choice, made throughout this paper, is a 
Dirac delta function so that T(u) = L 

iv) The total range of search for solutions that match the 
input tolerances on leakage or noise: K^in < Ka < 

^max ■ 

v) The maximum allowed number of interval bisections 
nbis of this K range in the search for solutions. 

There are additional, low-level, implementation-specific pa- 
rameter choices that affect the degree of numerical approx- 
imation made when calculating Aaij and Bai- These are 
described in Section 14.21 and the values used are tabulated 
(along with fixed choices of the soft inputs) in Table [T] 

The algorithm also requires input values and functions 
which are not chosen, and are determined by the images in 
question. These "hard inputs" are assumed to be fully known 
in advance, and are simply: 

i) The input images arranged into a continuous vector Ii 
and corresponding input pixel center positions r^. 

ii) The PSF at the position of each pixel center Gi (r). 

iii) The noise co variance of the image pixels Nij . 

Once all of these inputs are given the problem is fully spec- 
ified, and the algorithm can begin constructing the optimal 
linear transformation Tai- 

3.2. System Matrices from the Fast Fourier Transform 

The first significant computational task in the algorithm is 
the calculation of the system matrices Aaij and Bai, and the 
scalar Ca, given in equations (fT7])-(fT9l). We assume that Gi (r) 
and r(r) will be submitted to the software as fully-sampled, 
discrete images, although if some sufficiently accurate ana- 
lytic form is known this could also be used instead. 

A numerical approximation to the integrals of equations 
(fT7]i-(fT9]l might then be calculated directly for each Q!,i,j 
element in turn. This approach begins by constructing dis- 
crete representations of G'i(u) a nd r(u ) using the Discrete 
Fourier transform (DFT: see, e.g.. lPresse t al. 1992) of the in- 
put Gi(r) and r(r), then using these functions to take a di- 
rect sum over the integrand. However, a DFT of finite size 
in truth represents each PSF image as a periodic function in 
real space, which can causes erroneous successive maxima in 
the system matrices wherever the combination (r^ — r^) ex- 
ceeds the dimensions of the original PSF image; the same is 
true for (r,; — Rq) when calculating Bai- This problem can 
be overcome by zero-padding Gi(r) and r(r) to match the 
spatial extent of the input r,; and output Rq, and then taking 
the D FT. The Fast Fourier Transform (EFT: e.g., iPress et aP 
119921) algorithm performs this task highly efficiently. 

Yet this schema becomes prohibitively slow in the case 
of large, well-sampled PSF images. If we assume that the 
PSF images are better sampled than our input /; by a fac- 
tor npsF along each dimension, then the EFT of the zero- 
padded PSF images requires O [rip^pnlog (^.p^pn)] calcula- 
tions. That cost may be borne, but going on to calculate Aaij 
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then requires a further O (^psp^i^) calculations in total after 
each element is calculated in turn. This cost is high: for a 
single, small, 100 x 100 input image we have n = 10**, and 
if we assume that npsp = 3 we arrive at an estimated ~ 10^^ 
floating point operations required for the calculation of Agjj. 

Another approach is possible. Inspection of equations (fTTi 
& (fTsT l shows that we may write expressions for Adj and 
Bai directly in terms of inverse Fourier transform operations, 
as follows: 



G:(u)G,(u) (r,-r,), 



and 



= -2^-1 [G,(u)f*(u)] (r, - R„) 



(21) 



(22) 



The determination of A^ij and B^i thus proceeds by first 
constructing high-resolution lookup tables for these quanti- 
ties, using the inverse DFT of the G*Gj and G^r* arrays. As 
additional zero-padding may be added to these complex func- 
tions in Fourier space, the lookup tables may be sampled at 
essentially arbitrary spatial resolution (physical memory con- 
straints notwithstanding) following the inverse DFT back into 
real space. The value of each element in Aaij and Bai can 
then be determined b y high-order polynomial interpolation 
(e.g. iPress et al.|[T992 ) of the lookup tables to the locations 
Tj — Ti and Ti — Ha, respectively. 

Assuming once more that the PSF images are zero-padded 
to cover the spatial extent of the input images, the lookup 
tables will each require O [rip^pn log (ripgpn)] calculations 
using the FFT. A polynomial interpolation in two dimen- 
sions across regularly-gridded data is an operation of com- 
plexity 0(2npj,,y), where ripoiy is the degree of the interpolat- 
ing polynomial. Therefore, the determination of all elements 
of Aaij and Bai via lookup table requires 0{2n. 
0{2n. 



poly 



and 



poly 



calculations, respectivelyQ This process will 
take longer than the time required to construct the lookup ta- 
bles themselves, but the order n reduction in complexity seen 
over the direct-summation approach is a significant improve- 
ment. It is this technique that is used in the IMCOM prototype; 
it is noted that both approaches are trivially parallelizable, 
which is typically the case for calculations made throughout 
the method. 

3.3. Eigendecomposition of the Aq Matrix 

After the calculation of the system matrices of equations 
(fT7]i-(fT9]l. the next computationally challenging task is the so- 
lution of equation (l20l i to find Tai- This equation may be 
rewritten as 

Ta{Aa+K-N) = -^Ba (23) 

in matrix form; in the case of a single, fixed k for all a, 
this system is simply an m-element ensemble of transposed 
Mx = b matrix- vector equations. Here Alij = Aaij + uNij 
will be real, symmetric and positive definite, and each solu- 
tion (corresponding to a row of Tai) can be found efficiently 
and in parallel using the Cholesky decomposition M = LL"^ 
( iPress et al.lll992l:lAnderson et al.| [l999). This decomposition 
need be performed only once for constant k; subsequently, 

^ The authors note that the calculations presented in this paper were in fact 
made using a polynomial interpolation schema that allows arbitrarily-spaced 
points along the abscissa, but at a greater relative cost of 0{n'^^^^) operations 
per interpolation. 



only the back-substitution step need then be repeated m times, 
once for each row of ~Bai/2. The full solution for constant 



K therefore requires 0{rr 



n) ca lculations. 



However, as stated in Section \3A\ the goal of this imple- 
mentation is to allow the user to specify S™" or t/"^" for 
the output image, varying k from pixel to pixel to satisfy this 
requirement. If we assume that iteratively finding such a solu- 
tion requires rttdes attempts, then simply solving for Tai using 
Cholesky decomposition comes at a cost of 0(rttries'T^"^"^) cal- 
culations. We must seek to avoid fourth powers of n and to, as 
even for small images each of these numbers quickly exceeds 
- 10^ 

The problem may also be tackled by considering the eigen- 
decomposition of the system matrix Aaij, which is symmetric 
positive definite and may thus be written as 



Aa = QAQ^ = QAQ^ 



(24) 



where Qij is the n x n matrix whose columns are the 
n eigenvectors of Aaij, and the diagonal matrix Ay = 
diag ( Ai , A2, . . . , A„) contains the corresponding eigenvalues 
(e.g. ,Arfken & Weber.2005.) . The inverse is then given simply 
by 

A-i - QA-iQ^. (25) 

This fact can be utilized in the solution for Taij, if we as- 
sume that the input noise covariance is diagonal Nij = 
diag {{rjl), (77I), . . . , {Vn))- Whether this is a good or poor 
assumption will depend on the detectors used; off-diagonal 
terms mig ht occur due to co rrections for inter-pixel capaci- 
tance (e.g. Barron et al. 2007). It should be noted that even if 
Nij is not purely diagonal no systematic error is introduced 
in the output image, since the calculation of the leakage ob- 
jective Ua is unaffected. The method uses the input noise co- 
variance solely to minimize the output noise, and so imperfect 
assumptions about Nij will simply cause the resulting T^aa to 
be less than optimal. Any general Nij can later be propagated 
through to give the accurate output covariance matrix us- 
ing equation ^ and the solution for Tai as determined using 
the diagonal Nij approximation. 

Assuming a diagonal input noise Nij, we may then de- 
mand that the input image li is pre-multiplied by a scaling 

matrix diag {l/iviy^^ 1/(77!)'/', ■ • ■ , ^/ivlV^^) such that 
Nij = In, the n X n identity matrix. We then exploit a useful 
invariance property of the eigenvectors of a symmetric ma- 
trix under linear combination with the identity matrix: given 
equation ( l24b . it is simple to show that 



Aa + Kain = Q (A + 



(26) 



For our scaled image we may thus write 

(A„ + K„N)"' = Q(A + kJ„)"'Q^, (27) 

where the central matrix is purely diagonal and trivially cal- 
culated. 

The number of operations required for the initial decom- 
position of equation (l24l) is 0{n^), and this decomposition 
need only be perform ed once and may be parallelized (e.g. 
lAnderson et al.l [r999h . We show in the following Section 
13.41 that, subsequently, trial values of 'Saaina) or Ua{na) 
may be calculated for each output pixel a at a cost of 0{n) 
operations per trial solution. This implies a total cost of 
O [n^ + 2n^m + riingsnm) for calculating Taij, including an 
extra 0{n?m) overhead for necessary matrix manipulation as 
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will be discussed in the following Section. This represents 
a significant economy over the repeated use of Cholesky de- 
compositions and justifies the use of the eigendecomposition. 
This procedure also yields a figure for the condition number 
for the matrix A^ij, which gives a useful indication of the 
accuracy and stability of results regarding the input system. 

3.4. Solution for k by Interval Bisection 

Having calculated the eigenvectors Qij and eigenvalues 

that 



of the system matrix Aaij, we see from equation 
these may by used to trivially express the matrix inverse re- 
quired tocalculate in equation (l20l i. In AppendixlAl equa- 
tions (IA3b & (IA4b also give expressions for Sqq, and Ua in 
terms of this matrix inverse. Substituting equation (l27t and 
expanding, we find 



4^ 



. ^aj Qj 



and 



BajQji 



(28) 



(29) 



where here we have chosen to write the matrix multiplication 
explicitly for clarity. These expressions may be used to gen- 
erate trial values of Sqq (kq) or Ua {ko) directly, without first 
calculating Taij ■ The object 



(30) 



will be referred to as the projection matrix, and gives the 
scalar product of each row-vector Bq, with the successive 
eigenvectors Qij of Aaij ■ This projection matrix can be cal- 
culated in advance, immediately following the eigendecom- 
position of Aaij, at a cost of 0{ii?m) operations. Once this 
is done, each trial calculation of Y^aa or Ua made using equa- 
tions (l28T l & (|29] | comes at the low cost of 0{n) operations. 
After each is found, we may then use the projection matrix 
once more to write the final solution for T„, as 



T --^-^ 



-Paj Qi j 



2 ^ A, 



(31) 



at a further cost of 0{n^m) operations. 

The problem is then, for each output pixel a, to find the 
value of Ka that satisfies S™" - AE™" < < S™'' or 
jrjmax _ AC/™" < Ua < ^/™^ The simplest and most ro- 
bust algorithm for finding such a solution is interval bisection 
(|Press et al. 1992), and this is the method we adopt here, iter- 
atively bisecting some search range [n^in, i^max] at intervals of 
log K. If a root lies within this range the algorithm will find 
it, or will get as close as allowed by t he u ser-input maximum 
number of bisections ribis (see Section lTTl and Section l43] for 
some example choices of these parameters). 

Since there is some known directionality in the problem 
at hand, the algorithm may be adapted somewhat. In Ap- 
pendix |A] we show that the derivatives of J^aa or Ua with 
respect to k are everywhere negative or positive, respectively. 
Taking the case of Ua, this means that the first trial value 
should be calculated using K^in- if this solution does not sat- 
isfy Ua < i/™" then there can be no improvement on ^min 



and this value is adopted as Ha- For Sq-q, we likewise start on 
Kmax before proceeding to bisect. An output image of the in- 
terval bisection-derived Ka is also generated, so the user may 
see if this range needs to be expanded. In addition, it is also 
relatively fa st to generate an image of the Ua derivative given 
in equation ( IA6l l, which allows the user to judge convergence 
upon the fundamentally limiting lower bound, lim^^o+ Ua- 

The total cost of solving for Taij is therefore 0{n^ + 
2n^m + riaie^nm) as given in Section 13.31 which includes 
the eigendecomposition of Aaij, the calculation of the projec- 
tion matrix, the ntries trial solutions to Ka for each of m out- 
put pixels, and the final determination of Taij using equation 
yTJ- It can be seen that the dominant contributions come from 
the eigendecomposition of Aaij and calculating Pai, which 
lessens the need to seek more efficient (and potentially less 
robust) algorithms for the subsequent root-finding. 

4. AN ILLUSTRATIVE EXAMPLE 

We now demonstrate the use of the method by means of a 
worked example. As a primary goal of this work is the anal- 
ysis of survey images for a dark energy mission, we perform 
a simple simulation of observations made using an example 
design concept under consideration for the WFIRST mission. 
This particular design employs a D = 1.3 m primary mirror 
with an off-axis configuration for the secondary mirrors and 
detector array. This provides an unobstructed aperture, both 
increasing the flux throughput of the instrument and reducing 
the integrated light in the outer wings of the PSpQ 

We will simulate monochromatic observations at the near- 
infrared wavelength A = 1 /im, close to the center of the range 
proposed for the WFIRST imaging survey (Gehrels 2010). It 
should be noted that observations may likely take place at 
longer wavelengths than this value, and over a broad filter 
band. This will have two desirable effects, smearing together 
diffraction rings in the PSF and introducing a cutoff at lower 
spatial frequencies in the MTF. The A = 1/im monochromatic 
system adopted for these examples is, in this sense, more de- 
manding than we expect for a realistic dark energy mission. 

In Figure[T](left panel) we show a simulated PSF, generated 
using the ZEMAX software (Ed Cheng, John Lehan, David 
Content, & the WFIRST Project Office, priv. comm.), for this 
WFIRST design concept at this wavelength. The model in- 
cludes additional astigmatism to mimic aberrations due to fab- 
rication errors, and coma to simulate possible misalignments 
due to thermal drift and realistic imperfections in installation 
and testing. For the camera we assume an array of 18 /im 
Teledyne Hawaii-2RG (H2RG) detectors, corresponding to 
an image sampling of 0.18 arcsec. For the purposes of this 
demonstration we model the pixel response as a simple box- 
car function. 

We also mod el the effects of charge diffusion using re- 
sults from B arron et al.l (12007b . who measure a projected, one- 
dimensional diffusion length of Iqd — 1.87/im for a hyper- 
bolic secant diffusion function JcD(Aa;) cx sech(Aa;/ZcD), 
where Ax is the distance of the collected charge from the lo- 
cation of the electron-hole pair. Since we seek a de-projected 
charge diffusion model, some assumptions are unavoidable: 
for simplicity, we model the two-dimensional charge dif- 
fusion as a zero mean, circular Gaussian of width ctcd = 
7r/cD/2, preserving the variance (and hence diffusion length) 

* More details regarding tliis family of designs can be 
found in the report of tlie Interim Scien ce Working Group, 
http://wfirst.gsfc.nasa.gov/science/ISWG_R eport.pdfl 
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Figure 1. Left panel: Simulated PSF for a WFIRST mission concept with a 1.3 m unobstmcted primary at a wavelength of 1 fim. Right panel: Image of Gi(r) 
for this PSF, incorporating additional charge diffusion (CD) and the pixel response (PR) for 0. 18 arcsec detectors. Both images use a logaiithmic scale. 
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Figure 2. Left panel: Magnitude of the Modulation Transfer Function (MTF) for the WFIRST concept PSF of Figure [T] Right panel: Image of |Gi(u)|, the 
magnitude of the MTF corresponding to Gi (r); the sine modulation of the pixel response can be clearly seen at 1/0.18 = 5| cycles arcsec"^. Both images use 
a logarithmic scale. 



of the function. 

Convolving the WFIRST concept PSF with the pixel re- 
sponse boxcar and the charge diffusion function, we arrive 
at the functional form of Gi (r) for this demonstration, shown 
in Figure [T] (right panel). For this first demonstration, G,;(r) 
will remain constant throughout the input image. 

The Modulation Transfer Function (MTF) is defined as the 
Fourier transform of the telescope PSF. The function G'i(u) 
therefore represents the MTF conjugate to G,;(r), and the 



magnitude of this complex object can be seen in Figure |2] 
for the two corresponding PSFs of Figure [1] As can be seen 
from Figure ID the system is bandlimited at the fundamental 
frequency corresponding to D/X ~ 6.3026 cycles arcsec^^. 
Therefore, a sampling interval of \/2D = 0.079333 arcsec in 
the output image Ha is the require ment for c ritical sampling 
according to the sampling theorem (Marks 2009). This output 
sampling rate is therefore of fundamental interest in demon- 
strating that survey design and multiple dithering strategies 
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will allow a dark energy mission to overcome an undersam- 
pling detector array in the focal plane, and thereby utilize all 
available spatial information for weak lensing measurements. 
Hence, X/2D = 0.079333 arcsec will be adopted as the de- 
sired output sampling interval for all the demonstrations in 
this study. 

4. 1 . Input Images: Random Dithers, Random Rotations 

From the WFIRST concept described in Section |4] it has 
been possible to create a model for Gi{r) that shares many 
of the properties expected for a future dark energy mission. 
However, this represents only one of the input quantities de- 
scribed in Section im not least, input images themselves are 
required. 

As a demonstration of the technique's flexibility we will at- 
tempt to reconstruct a fully-sampled output image using a set 
of three, randomly dithered and rotated images. The spatial 
configuration r, of the three dithered exposures can be seen in 
Figure [3] each of dimension 20 x 20 pixels^. Also shown are 
the corresponding sample locations of the desired, fully- 
sampled output image. It should be stated from the outset that 
these three dither patterns are not expected to be sufficient to 
achieve full sampling, but are being used deliberately to show 
the behavior of the algorithm when unable to recover certain 
output pixels. The user may of course use images with any 
pixel locations, there being no requirements of regularity or 
homogeneity in the formalism. 

For the physical image on the sky /(r) we choose a sin- 
gle galaxy, modeled as a simple elliptical, exponential pro- 
file. The profile is given an axis ratio a/b = 3/\/5, a semi- 
major axis scale length (aligned along the x direction) of 
a = 3/ (4\/5) arcsec, and a central peak surface brightness 
of 50 (in arbitrary units). This model galaxy is shown in Fig- 
ure |4]in high-resolution. 

This image is then convolved with the GJr) of Figure [1] 
and sampled at the locations shown in Figure|3]to generate the 
array of input images that make up I^. Since each image is ro- 
tated, before calculating Aaij and Bai the PSF Gi (r) for each 
exposure must be rotated into a co mmon coordinate system in 
which r(r) is defined (see Section l4~2l below). In this exam- 
ple we choose r(r) to be simply the unrotated Gi (r) of Figure 
[H so that only the input PSFs need be rotated before solving 
the system. For precise comparison of output and input im- 
ages we do not add stochastic noise to li, but will nonetheless 
assume a unit diagonal noise covariance Nij = I„ to illus- 
trate the behavior of Y,aa in a more realistic, noisy imaging 
scenario. 

With Gi (r), li, Tj an d Nij we have provided the three "hard 
inputs" of Section im but there remain choices to be made re- 
garding the values of user-specified input quantities. We adopt 
a simple Dirac delta function for T(r) so that T(u) = 1. 
This fixed parameter choice is given in Table[Tl where all such 
implementation-specific parameter values are given. The re- 
maining choices that must be made relate to the process of 
finding solutions for L/"™ or S™", and so we defer a discus- 
sion of these input quantities until Section [43] We first illus- 
trate the system matrices A^ij, Bai, which may be calculated 
immediately. 

4.2. System Matrices 

Taking the input functional forms of Gi{r) and r(r), and 
the pixel locations and R^, we may calculate Aaij and B gj 
using the FFT-based prescription described in Section U!2l As 



the accuracy of these objects is crucial, and their computation 
costs are relatively low in comparison to later operations, it is 
desirable to generate high-resolution lookup tables for Aaij 
and Bai and then interpol ate u sing a high-order polynomial. 

As discussed in Section |372l the images of Gi(r) and r(r) 
are first zero-padded so that their spatial dimensions (in phys- 
ical units) equal or exceed the maximum spatial extent of the 
input images that make up li. It is noted that the WFIRST 
design concept PSF of Figure [T] is supplied at relatively high 
resolution (npsF — 11: see Section 13.2b . For rotated input 
frames, each Gi (r) must be rotated accordingly to the com- 
mon coordinate system in which r(r) is defined. Each image 
Gi{r) is rotated using Lagrange polynomial interpolation at 
order Upoiy — 7. Then Gi(u) and r(r) are calculated using 
the FFT algorithm. Limited only by the memory available 
(2GB in a 32-bit architecture), we may then add further zero- 
padding to the bandlimited Gi{u) and r(r). We choose to 
zero-pad each such that their linear dimension increases by a 
further factor ripad compared to the starting size of Gi{u) and 

r(r) (which were already padded once in real space, up to the 
size of the input images). In this example, and throughout this 
Paper, we choose ripad = 3 (see Table[T]i. 

The inverse FFT is then used to build the lookup tables de- 
scribed by equations (l2Tl i & (|22] | which, following the zero- 
padding, have npsF x ripad entries along the length of each 
input pixel in each dimension. Finally, these tables are in- 
terpolated using a polynomial at order npoiy = 7 to estimate 



values of Aaij and Ba 



at the locations of the vectors r, — r,; 



and Ti — Rq, respectively. The resultant matrices can be seen 
in Figures 12] & [6] 

Once Aaij and Bai are calculated, the final stages in the 
preparation for the solution of Tai are the eigendecomposi- 
tion of the Aaij matrix to give Xi and Qij, and the subse- 
quent calculation of the projection matrix Pai- In Figures |7]& 
[8] we plot the eigenvalues and lowest-order eigenvectors for 
Agii in this exarnple. T he condition number H (see, e.g., 
IChenev & KincaidllI998l) for a symmetric, positive-definite 
matrix such as Ar_ 
values as 



iQ,jj is simply expressed in terms of the eigen- 

max[Ai] 



^ {Aaij) — 



min[Ai 



(32) 



A system is said to be singular if the condition number is infi- 
nite, and typically ill-conditioned if log^o (^) ^ the numerical 
precision of at which each matrix element is stored. For the 
illustrative example in this Section, the condition number is 
found to be H (Aaij) — 8.369 x 10^, which implies the sys- 
tem is adequately well-conditioned if using double precision 
arithmetic with machine epsilon e = 2^^"^ = 1.11 x 10~^^. 
In practice the matrix that must be inverted is Aaij + HaNij 
so, given appropriate choices of Ka, the method may proceed 
even when H is large. Nonetheless, calculating H in this way 
provides both a useful check on the eigenvalues (which 
must all be positive since Aaij is a symmetric positive def- 
inite matrix), and an idea of when solutions for Tai might 
become unstable for small Ka ■ 

The projection matrix Pai is then constructed from the 
eigenvectors Qij and the matrix Bai via matrix multiplica- 
tion as shown in equation ( l30l l. At this stage the problem is 
then fully set up, and we are ready to begin the solution for 

4.3. Output Images 
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Figure 3. Left panel: Input pixel locations for three randomly dithered and rotated exposures of the same camera. Red, green and blue crosses represent each 
of the three exposures. Right panel: Output image sampling locations (black crosses), placed at a spatial interval of 0.079333 arcsec so as to fully sample 
the WFIRST concept image plane. 



Table 1 

Fixed, IMCOM implementation-specific parameters used in calculating the results of Sections |4]&ll] 



Parameter 


Value 


Description 


T{r) 


S{r) 


Kernel used for calculating leakage objective Ua (see equation[TT). 


"pad 


3 


Multiple by which Gi (u) and r(u) sue further zero-padded in u when calculating 




Acij and Bai (Section[4.2r 




7 


Polynomial order used when rotating Gi{r), and interpolating Acdj and B^i from 


l.llCa X 10-16 


lookup tables (Section[4.H. 


'^min 


Minimum Ka for interval bisection search, set by machine precision 
(Sections [3l&|43l. 




9.01Ca X lO^s 




Maximum k„ for interval bisection search, set by machine precision 
(Sections[3l&l43t. 




53 


Maximum number of interval bisections allowed in solving for each Ka 




(Sections [3l&l43t. 





10 


20 


30 


40 



3.0 - 
2.5 - 




1.0 - 
0.5 - 

0.0 r ' 

0.0 0.5 1.0 1.5 2.0 2.5 3.0 
X [arcsec] 

Figure 4. Surface brightness /(r) used in Sections|4]&|5] a simple, idealized 
model of an elliptical galaxy. 

T he so lution of equation (l23T l for Tai as described in Sec- 
tion |33] proceeds in two separate stages. First, values of Ka 



for each output pixel are independently determined so as to 
attempt to keep Ua or Y,aa within the pre-set limits. Second, 
and at far greater relative computational cost, these values are 
substituted into equation ( |3T] i to calculate Tai- 

To determine Ka, it is necessary to complete the set of 
choices described in Section 13.11 For the purposes of this 
example, we will demonstrate the solution of Ka both by 
setting requirements upon Yiaa and by setting requirements 
on Ua- For both cases, we search within the broad range 
r^min^^max] ^ [lAlCa X 10-l^ 9.01(7„ X 10^^] (scc Table 
U^, effectively being limited only by the machine precision e. 
We also allow a total of nbi.s = 53 interval bisections to find 
Ka- Computational costs are negligibly affected by choices 
for these parameters. 

In Figure |9] we plot results from the iMCOM code having 
required E™" - AE™" < Eaa < S™' ' where E™^" = 1 and 
AE™" = 10^^. As the input noise covariance is in this case 
specified as a unitary, diagonal Nij = I„, this corresponds to 
requiring that the output image Ha maintain the same level 
of noise variance as compared to the input image pixels. Fig- 
ure shows the value of Ka selected at each output position 
Rq in order to fulfill these requirements on E„q, and Figure 
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Figur e 5. System matrix Adj for the image configuration of Section|4]and 
Figure [To] 
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Figure 8. The first tliree eigenvectors of the A^ij matrix shown in Figure|5] 




Figur e 6. System matrix Bdi for the image configuration of Section |4]and 
Figure [Tol 
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Figure 7. Eigenvalues of the A^ij matrix shown in Figure [3] The matrix in 
this example has a condition number of K = 8.369 X 10*^ . 



|9}3 shows the variation in the resultant leakage objective Ua 
across the output image. This latter is plotted in it s natural 
units of Ca (see Appendix lAl equations I A8 1 & |A9l ) . Degra- 
dation in image fidelity, manifesting as increased Ua, can be 
clearly seen in the corner regions of the output in Figure 0) 
where input pixel coverage is reduced, but also within central 
locations where chance alignments of the input pixel pattern 
lead to encircled regions of sparse sampling. 

Overall we see that, given the requirements on the noise 
properties, the randomly dithered, three exposure system pro- 
duces output with Ua ^ O.OlCa or worse over a significant 
fraction of Ha- This is not the level of fidelity that will be 
required for a dark energy mission, but the demonstration il- 
lustrates how aspects of the input dither pattern impact the 
quality of the recovered output. 

The corresponding noise variance map Y,aa can be seen in 
Figure |9};. As shown in Appendix [A] the output Eqq, can 
be made arbitrarily small by increasing Ka, but not necessar- 
ily arbitrarily large as Ka decreases. The map of Y^aa shows 
where occasional chance alignments of multiple input pixels 
produce noise that is lower than S™™ — AS™", and compar- 
ison with Figure shows that these points indeed lie where 

«• — „min 
l%a — n-Q ■ 

In Figure |9}l we compare the iMCOM-reconstmcted image 
Ha to an independently-generated target image J^, calculated 
using equation ^ and the highly-oversampled /(r) as shown 
in Figure|4](the code for generating Ja is in fact written in the 
IDL language and shares no routines with iMCOM). These 
two image vectors are used to calculate the squared fractional 
residual, which we define as 



j-,2 {Ha J a) 

= 



(33) 



Values of depend upon /(r) directly and so cannot them- 
selves be used to locate an optimal Tai, but in simulated 
examples they provide a useful comparison: for a well- 
controlled system we typically expect < Ua/Ca [al- 
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X [arcsec] X [arcsec] 



Figure 9. Map of solution (a), output leakage objective Ua (b), noise variance Sea (c), and corresponding (d) for the input choices S™^ = 10 ^, 
AS™^ = 10~^ in the demonstration of Section|4] The color scale in the map for Ua does not reflect the full range of these values, but is selected to most 
clearly illustrate the dependence of Ua on the input dither pattern r ^ . 



though this may be violated in practice if Ja ~ or /(r) con- 
tains large or abrupt intensity variations such as those caused 
by bright stars or cosmic ray impacts]. The variation in 
across the output image plane for this example can be seen in 
Figure |9}l- This shows that < Ua/Ca across much of the 
output in this example but can vary sharply due to the large 
variations in Ua, particularly towards the edge regions of the 
image where /(r) is smaller 
We now examine a case in which the desired output im- 



age properties are specified in terms of the leakage objective 
Ua rather than the noise variance. In Figure we plot re- 
sults from the iMCOM code having required J/™''" - AC/™^" < 
Ua < C/^^ where J/™" = IQ-^Ca and AJ/™" = IQ-^C^. 
It can be seen from Figure [TOb that in fact Ua > IO^^Cq 
in the sparsely-sampled corner regions of the output region. 
Here in the corners the system has saturated at the lower limit 
for Ua within the broad range [k™" , k™""] , the practical mini- 
mum possible using double precision arithmetic and possibly 
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Figure 10. Map of Ka solution (a), output leakage objective Ua (b), noise variance E^a (c), and corresponding (d) for the input choices {7™^" = 10 
At/™"" = 10~* in the demonstration of Section|4] The failure to reach Ua < 1 x 10~^Ca in comer regions, even for = k™", can be clearly seen; the 
localization of peaks in the noise variance around encircled zones of reduced pixel coverage is also apparent. 



approaching the theoretical minimum given by equation (|A9) . 

Figure [Tob shows the output Y,aa for this solution: rapid 
variation is seen in the noise properties for this output image, 
with factors of ^ 10-100 greater variance than the input Nij 
seen in some regions. These include the comer regions but 
also, as was the case for Ua in the previous example, encir- 
cled regions of sparse sampling that are reflected in the output 
Eqq. There are also similar traces of these effects in the 
map of Figure [Tol l. although these are most marked in the 



edge and corner regions where sparse input sampling com- 
bines with small values of Ja in the target image. 

Overall, it is seen that restricting Ua < IO^^Cq, (a rather 
modest ambition if attempting precision scientific work) is 
difficult for the input image configuration and PSF of this ex- 
ample system. While, as shown in Figure |9] it was possible to 
create an image with benign noise properties, this comes at a 
cost of significant Ua across much of the image. Large values 
for Ua require significant leakage LaCRa — r), and therefore 
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imply a non-negligible, biasing residual (Za). The input sam- 
pling is too sparse to allow stable, unbiased reconstruction at 
this output resolution. While this would have been apparent 
from the outset to someone with experience in stacking imag- 
ing data, the linear formalism allows it to be quantitatively 
diagnosed using and Ua in a way that is independent of 
/(r). This capability makes the method useful as a survey 
design tool. 

In this Section we have attempted to demonstrate the flex- 
ibility and power of the linear image combination formalism 
developed in Section |2] by tackling a difficult case. For the 
three randomly-dithered exposure example, where the input 
sampling turns out to be insufficient for a fully-sampled, un- 
biased output, we have illustrated how the objective quantities 
Ua and are used to understand the properties of Ha- In 
placing simple, user-specified constraints on these objective 
quantities we have shown how the iMCOM implementation of 
Section|3]can be used to efficiently explore the optimal trade- 
off between noise and fidelity, using a freely-varying Kq. In 
trickier cases such as that presented here, there are further 
synthetic options that might be used to improve noise and bias 
in Ha'- specifyi ng r (r) as a broader PSF is one possibility 
(see also Section l53] l. 

For the purposes of designing a dark energy mission, how- 
ever, it is useful to address the problem in a slightly different 
manner. Given a telescope PSF and an undersampling focal 
plane detector configuration, we may ask how many dithered 
exposures are needed, and in which pattern, to generate an 
unbiased, fully-sampled output Ha- This represents an im- 
portant practical application of the techniques developed in 
Sections |2] & |3] and it is using this capacity as a design tool 
that we explore a range of dither patterns in the following Sec- 
tion. 

5. TESTS OF MULTIPLE OBSERVING SCENARIOS 

In this Section we assess the ability of a number of dither 
pattern configurations to reconstruct fully-sampled images 
with a low level of leakage, a sure requirement for the pre- 
cision photometry and shape measurement that will be neces- 
sary for the imaging component of a dark energy mission. For 
the purposes of these tests we define low leakage as having 
Ua < lO^^Ca, implying control of the output PSF proper- 
ties to one part in 10"*. This ensures that uncontrolled changes 
to the PSF from linear image reconstruction will account for 
only a fraction of the total systematic err or budget for a dark 
energy sui-vey ( Am ara & Refregierl 1200 8 ). We then exam- 
ine the capability of the WFIRST design concept described in 
Section!?] with an optical PSF given by Figure [T] and a focal 
plane of H2RG detectors sampling at 0. 18 arcsec, to generate 
fully-sampled outputs that are unbiased at this level of leak- 
age. 

The problem is set up as follows: We adopt C/™^" = 
IO^^Cq and AL/™™ = 10"^"Cq as our requirements on the 
properties of the output image Ha, whereas other inputs such 
as the undersampling input pixel scale (0.18 arcsec), G,;(r), 
r(r), and the parameters of Table [T] are all as described in 
Section |4] unless explicitly stated as otherwise- Using a re- 
alistic WFIRST design concept in this way, we can begin to 
explore what survey design strategies may be needed to over- 
come undersampling for a dark energy mission of this sort. 
We continue to set Nij = I„ for the easy interpretation of 
Eq,q results. 

The output Rq, will again be placed at an interval of 
0.079333 arcsec, the requirement for critical sampling. In 



each following test we then vary the input pixel sampling in- 
formation, adopting different dither patterns, numbers of ex- 
posures or detector coverage: in terms of the formalism, these 
variations are manifested as changes in and li- In this way 
we seek to show how the linear image combination formalism 
can be used to inform survey design, and make recommenda- 
tions about survey dithering strategies for this WFIRST design 
concept. 

5.1. The 2x2 Ideal Dither 

In Figure[TT|we plot results for four 36 x 36 pixeP input ex- 
posures, configured in a 2 x 2 ideal dither pattern. The upper- 
left panel. Figure [TTk. illustrates this dithered configuration of 
input exposures (colored points) along with output sampling 
locations (black crosses), focused in upon a small sub-section 
of the total image area for clarity. Figure ITTb shows the mini- 
mum leakage objective Ua that may achieved at fully-sampled 
resolution for this dither pattern; as can be seen, the threshold 
value of J/™" = 10~®Cq is often exceeded (it is reached only 
where output pixels closely align with input pixel locations, 
indicated by the regularly spaced minima in [/„)■ 

The map in FigurefTTl: shows the output noise variance Y.aa 
for the 2x2 dithered image. Interestingly, this is both reason- 
ably stationary (i.e. approximately constant), and comparable 
in amplitude to the input noise variance, across a large area of 
the output. We note that in this Figure, and in many subse- 
quent plots of EoQ, we choose the color scale not to show the 
full dynamic range of J^aa but to better highlight behavior in 
the central regions of the image. 

At the extreme edges of the output, where input information 
is most scarce, may exceed 10^. The size of such noisy 
edge regions is primarily determined by the extent of the PSF 
Gi (r) within each contributing exposure. However, as will be 
discussed in Section|6] these regions are not likely to be used 
in a realistic survey strategy: Figure (Tsjillustrates that it is not 
necessary for these edge regions to fall within the output Ha 
at all. Instead, by tessellating successive Ha that cover only 
the central regions of overlapping input patches, a contiguous 
output image can be constructed for which all pixels are suffi- 
ciently far from an input edge. In this Section, however, these 
noisier edge regions are kept so as to demonstrate the effect. 

Figure [Till shows the map of across the image. The 
ciTiciform pattern seen in F^ is due to the light profile of the 
input /(r) seen in Figure |4l this fact was verified by varying 
the centroid position of /(r), which was seen lead to equiva- 
lent changes in the position of the cruciform. It can be seen 
that F^ < Ua/Ca across much of the output for this semi- 
realistic galaxy image. 

5.2. The \/E X \/E Ideal Dither 
In Figure [12] we plot results for five 32 x 32 pixel'^ input 
exposures, configured in a \/5 x \/5 ideal dither pattern. The 
upper-left panel. Figure [T2h . illustrates this configuration of 
input exposures for a small sub-section of the total image 
area. This is an interesting dither pattern which arranges five 
exposures per unit cell in a slanted, but regular, grid config- 
uration. The pattern has chirality: a reflection about the line 
y = X produces a different, but functionally equivalent, set of 
dithers. 

Figure [T2b shows the Ua achieved for this -\/5 x ^/E dither 
pattern: across the broad central region of the output image 
the requirement Ua < lO^^Ca is now successfully met. Fig- 
ures [T2b & d show the noise variance Eqq and squared frac- 
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Figure 11. Inpu t pixel configuration (a) and leakage objective Ua (b), noise variance S^q (c), and squared fractional residual (d) for the 2x2 ideal dither 
(see Section lsTt . The color range chosen for Saa does not reflect the full range of values in the output, but highlights variations in the image center 



tional residual F^. The former is seen to be approximately 
stationary, and comparable to the input noise variance, across 
the central regions of the output. The values of again form 
a vaguely-cruciform pattern that moves to align with the cen- 
troid of /(r) when this is shifted. Overall, the x i/5 dither 
pattern is able to successfully reconstruct a fully-sampled out- 
put image for this WFIRST design concept while preserving 
the input convolution kernel Gi{r) at a level Ua < lO^^Ca- 
A question remains as to whether a pattern such as this can be 
used to produce fully-sampled data in the presence of image 



pixel losses due to, e.g., cosmic rays, hot pixels, and other de- 
fects, or in the presence of variations in th e foc al pla ne plate 
scale: these issues are explored in Sections [5. 4I & I53] 

For illustrative purposes, we show in Figure [13] an exam- 
ple of an input exposure and the iMCOM-generated output 
for the -y/5 x \/5 ideal dither, where actual stochastic noise 
rji was added to the input pixels li before processing. The 
noise added was Gaussian with unit variance. The output Ha 
clearly shows the effects of increased J^aa in edge regions, 
but demonstrates the successful image reconstruction at the 
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Figure 12. Inpu t pixel configuration (a), leakage objective Ua (b), noise variance S^a (c) and squared fractional residual (d) for the \/5 X \/5 ideal dither 
(see Section lS^ . The color range chosen for Sqq does not reflect the full range of values in the output, but highlights variations in the image center 



image center. 

5.3. The 6 Exposure Random Dither 

In Figure [14] we plot results for six 32 x 32 pixel^ input 
exposures, configured in a randomly offset dither pattern; un- 
like in Section IJthere is no random roll in this example. Fig- 
ure [T4h illustrates this configuration of input exposures for a 
small sub-section of the total image area, along with the fully- 
sampling output locations. The 5 exposure random dither was 
also tested, but produced poor results with Ua > IO^^Cq, 



across much of the output in a small suite of randomly offset 
trial configurations. One such test, not atypical of the small 
set of trials, produced J^aa > 10 for 65% of the total output 
and > 100 for 22% of the total output. 

However, the leakage objective for the 6 exposure random 
dither in Figure [T4b shows good agreement with stated re- 
quirements across a large fraction of the output sample loca- 
tions. The Yiaa variance in Figure [l4b is also now signifi- 
cantly smaller than that which was seen in the five exposure 
dithers: an impressive impact from only a single additional 
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Figure 13. Left panel: One of five noisy input images dithered using the X \/5 ideal dither pattern of Section l531 constructed by convolving /(r) with d (r) 
and adding a stochastic noise rn to each pixel li. The noise added was drawn from a Gaussian distribution with unit variance and zero mean. Right panel: The 
combined output image Ha ■ Both images show the logarithm of the absolute value of each li and Ha , to most clearly depict the central, edge, and background 
noise-dominated regions of each image. 



exposure. As was done for the 5 exposure dithers, these re- 
sults were repeated internally for a small set of random 6 ex- 
posure dither patterns, and the results shown are typical of 
the results achieved. It should be noted that the noise Y,aa is 
somewhat non-stationary, varying with Rq, across the output. 

As has been seen, good results have been achieved for 
a controlled reconstruction of fully-sampled images using 
certain configurations of input WFIRST-like images. The 
-\/5 X y/E ideal dither performed well, as did a single realiza- 
tion of a random dither of six exposures (at the cost of slightly 
inferior noise properties). While the results shown for the ran- 
dom patterns were typical of a small internal test set of such 
patterns, the repeatability and range of these random dither 
results should be explored in more detail in subsequent work. 
Another investigation, which is within the scope of this Paper, 
is the robustness of linear image combination under identified 
image or pixel defects. This is now explored. 

5 .4. Cosmic Ray Impacts and Bad Pixels 

We perform a very simple test of linear image combination 
under the presence of known image defects, such as might be 
caused by cosmic rays, dead pixels, and similar phenomena. 
Tak ing t he input images for the \/5 x Vs ideal dither of Sec- 
tion |5.2| we randomly select 5% of input pixels to be flagged 
as not for use in the reconstruction of the input images. This 
equates simply to a contraction of the Aaij and Bai by re- 
moval of the afflicted pixel rows and columns. 

In Figure [T5b we plot the map of k.^ solution values along 
with, in Figure [T5b . the leakage objective Ua that is main- 
tained by the algorithm for this test. Also shown are the loca- 
tions of the 5% of input pixels that were thrown out, as small 
black points. Features in Kq, can be clearly seen where the IM- 
COM algorithm has attempted to satisfy requirements on Ua 
where input image information is scarce. 



The increased noise variance in Sqq, for these bad pixel 
locations can be seen in Figure [T5t. Indeed, the noise can 
be large in regions of missing input information, particularly 
where pixel removals are clustered. This information is rep- 
resented in the output covariance given by equation (|9]l, 
and is available to the user for weighting or fitting purposes 
if necessary. If missing pixels are expected to be clustered 
(such as around a single cosmic ray hit) the rapid increase in 
noise for such events may represent a cause for concern that 
should be investigated using simulations of greater realism. 
However, the level to which Ua (and indeed F^: Figure [TStl) 
may be controlled is perhaps more important for many appli- 
cations, and here results are encouraging. 

In Figure [16] we show equivalent results after randomly re- 
moving 5% of input pixels for the 6 exposure random dither 
discussed in Section l53] Once again, results for Ua are robust 
despite bad pixel losses, with Y^aa able to absorb uncertainty 
due to these small, localized instances of missing informa- 
tion. As for the \/5 x \/5 ideal dither, Sqq, may become large 
where missing pixels are clustered and this behavior should 
be investigated further 

It has been shown that the WFIRST design concept is able to 
produce output that is relatively robust to 5% pixel losses from 
random locations in the input images li, both for the \/5 x \/5 
ideal dither and a 6 exposure random dither configuration. 
Tolerances on Ua may be maintained at the cost of increased 
T^aa in regions of missing information: qualitatively this is 
the desired behavior However, there are other practical is sues 
that may be encountered for certain survey strategies, one of 
which is the possibility of significant plate scale variations 
between exposures when combining wide angle-dithered im- 
ages. We now investigate linear image combination in such 
circumstances. 
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Figure 14. Input pixel configuration (a), leakage objective Ua (b), noise variance (c) and squared fractional residual (d) for the 6 exposure random 
dither (see Section l531 . The color range chosen for T,aa does not reflect the full range of values in the output, but highlights variations in the image center. 



5.5. Focal Plane Plate Scale Variations 

Unavoidable geometric field distortions, due to optical 
aberration and the alignment of each detector array in the fo- 
cal plane assembly, are manifested as variations in the plate 
scale (angular pixel scale) across a telescope field of view. 
These geometric distortions may be determined using preci- 
sion astrometry, but it is not clear what impact they may have 
on a linear image reconstruction from multiple dithers. Typi- 
cally the distortion is smoothly varying, so that close pairs of 
pixels on the detector plane show smaller relative scale vari- 



ations on average than widely-separated pairs. With wide- 
angle slew dither strategies being considered for the WFIRST 
imaging survey, the behavior of linear image combination 
when facing plate scale variations in the input images is an 
important design consideration. 

We investigate the effect of plate scale variations using the 
6 exposure random dither pattern discussed in Section 15.31 
For each input exposure we vary the relative fractional plate 
scale by a random amount drawn from a Gaussian distribution 
of width 1%. The variation is applied to all pixels in each 
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Figure 15. Map of Kc solution (a), leakage objective Ua (b), nois e vari ance S^a (c), and squared fractional residual (d) for the \/5 X \/5 ideal dither after 
losing a randomly-selected 5% fraction of input pixels (See Section is!?! . Input pixels that were not used are indicated using small black crosses. The color range 
chosen for T,aa does not reflect the full range of values in the output, but highlights variations due to bad pixels in the image center 



given exposure, and is the same in both x and y directions: the 
effect is to uniformly expand or contract the pixel grid. The 
charge diffusion and pixel response component of each d (r) 
is appropriately modified to reflect the new input pixel scale, 
and r(r) is chosen as the Gi{r) corresponding to exposure 
with the largest pixels. 

In the upper panels of Figure [17] we plot the Ua and Eq,q 
that can be achieved for this dither configuration, having stip- 
ulated J/™" = IQ-^Ca and AC/™" = lO-^°Ca as usual. 



It is clear that the impact of ^ 1% plate scale variations is 
severe: Eq,q sho ws si gnificant increases as compared to the 
results of Section |53] whereas degradation Ua is more mod- 
est but perceptible. An increase in noise variance by a factor 
> 10 compared to the input pixels, as seen here over much 
of the output, is not a tolerable cost for fully-sampling a dark 
energy imaging survey. The success of linear image combi- 
nation is clearly sensitive to small variations in the input plate 
scale, and a mitigating strategy must be found. 
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Figure 16. Map of Kq solution (a), leakage objective Ua (b), noise vari ance Sqq (c), and squared fractional residual (d) for the 6 exposure random dither 
after losing a randomly-selected 5% fraction of input pixels (See Section l5!4l . Input pixels that were not used are indicated using small black crosses. The color 
range chosen for Saa does not reflect the full range of values in the output, but highlights variations due to bad pixels in the image center. 



One such mitigating strategy is to increase the number of 
dithered exposures used in the reconstruction. This was at- 
tempted, adding additional exposures to the 6 exposure ran- 
dom dither, up to a total of nine. Results were modest, show- 
ing only marginal improvements in I]q,q as extra dithers were 
added. 

Another option to mitigate the effects of scale variation is 
to exploit the user's freedom in setting the desired output PSF 
r(r). In the lower panels of Figure [17] we plot results us- 



ing the same input pixel configuration as in the upper pan- 
els, but having added an additional Gaussian smoothing with 
a — 0.09 arcsec (0.5 input pixels) to r(r). Dramatic improve- 
ments to the noise properties Sqq, can be seen. Ua is also 
improved enough to bring it below IO^^Cq for a far greater 
proportion of the output area (see Figure [TtI). 

We see, therefore, that the damaging impact of plate scale 
variations can be suppressed to some extent by the addition of 
synthetic smoothing to the desired PSF r(r). Considering the 




Figure 17. Upper panels: Leakage objective Ua (a), and noise variance Sea (b), for a 6 exposure random ditlier after adding ~1% variations in the pixel scale 
for each input exposure. Lower panels: Ua (c), and noise variance Saa (d), for the same input pixel configuration, but adding a Gaussian smoothing of cr = 0.09 
arcsec (0.5 pixel) to the desired output PSF r(r) (See Section|53). Plotted color scales do not necessarily reflect the full range of Ua or Scvcv, but are chosen to 
emphasize the central regions of greatest interest. 



effect upon the MTF r(u), this can be seen simply as a sup- 
pression of the high-frequency modes that are most challeng- 
ing to recover in the presence of plate scale variations. The 
choice of optimal smoothing filter will, in general, depend 
closely upon the nature of the plate scale variations encoun- 
tered between input exposures, and on the dither strategy em- 
ployed. The successful linear reconstruction of input images 
has been shown to be sensitive to this effect: further study 
aimed at finding optimal mitigating strategies, using more re- 



alistic simulations, will be needed. 

This concludes our first set of trials for the optimal linear 
image combination formalism. These early results suggest 
that the technique has merit as a design tool and can be used 
to explore competing dither strategies for generating over- 
sampled output. It has been used to identify encouraging ro- 
bustness of linear reconstruction to the presence of randomly- 
located, missing input pixels, but has also highlighted a sensi- 
tivity to plate scale variation that should be explored further 
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Figure 18. One approach for efficiently processing an entire dark energy 
survey is to split it into small patches as shown here. In this approach, in- 
put pixels li across multiple exposures are drawn from within a square of 
dimensions I X I (dashed line). From these, an output image Ha is recon- 
structed within a smaller, contained square of dimensions s X s (solid line). 
Successive patches of sky may be aligned so that the smaller output regions 
tessellate to cover the entire survey. 

Clearly, however, an ultimate goal for the method is also 
its use in constructing the oversampled images from actual 
all-sky imaging data supplied by a mission such as WFIRST. 
Therefore, we now turn to a discussion of whether such an ap- 
proach represents a computationally feasible (i.e. affordable) 
option for analyzing data from such a mission. 

6. PRACTICAL CONSIDERATIONS & 
COMPUTATIONAL COSTS 

The linear image combination formalism provides a natu- 
ral framework for constructing oversampled images, but an 
equally natural question begs itself: can the method be ap- 
plied to a 10 000 deg^ imaging survey (such as proposed for 
WFIRST) within an acceptable timescale, and using reason- 
able computing resources? To answer this question, we must 
begin by estimating the total number of floating point opera- 
tions required to apply the algorithm to a survey of this area. 

As discussed in Section |3] the solution for any given patch 
of sky requires 0{n'^ + 2n^m) floating point operations, ne- 
glecting terms that are lower order in n or rn. Therefore, the 
optimum data processing strategy is simply to split any imag- 
ing survey into as many small patches as possible: this is clear, 
since n for each patch is proportional to its area. Reduced 
memory costs and the ease of coarse-grained parallelization 
are other advantages of processing many small patches of sky. 
In Figure[T8]we illustrate one such patch schematically, noting 
the choice of different size input and output regions to avoid 
edge effects. The output regions (solid line) for such patches 
can in principle be laid down in a tessellating pattern to cover 
an entire survey area, with the larger input regions (dashed 
line) therefore overlapping to provide approximately uniform 
coverage free from the edge effects discussed in Section ISTI 

We will now consider the computational cost of such a pro- 
cedure. Using the labels of Figure [18] and taking the same 
WFIRST input pixel size and output sampling rate (0.18 and 
0.079333 arcsec, respectively) as Sections |4]&|5] let us con- 
sider patches of outer side I — 25 x 0.18 arcsec — 4.5 arcsec 
and inner side s = 38 x 0.079333 arcsec — 3.014 arcsec. 
This is a small patch, and thus an efficient means to cover 
a large survey, but is of sufficient size relative to the PSF to 
avoid significant degradation in Ha due to edge effects (see, 
e.g.. Section ISTTT l. A 10 000 deg^ survey requires 1.4 x 10^" 



such patches for full coverage. Assuming five dithered input 
exposures gives n — 3.125 x lO'^ and m = 1.444 x 10'^ for 
each patch, resulting in order + 2n^m ~ 5.8 x 10^° float- 
ing point operations per patch. To process all patches inde- 
pendently to cover the total survey area thus requires ^ 10^^ 
floating point operations, a considerable computational cost. 

However, because of the extreme parallelizability of the 
problem the use of supercomputing resources is a realistic 
option. The RoadRunnefl supercomputer at Los Alamos Na- 
tional Laboratory is capable of a sustained rate of 1 PFLOPS 
(10^'^ Floating Point Operations Per Second) in double pre- 
cision arithmetic, and would therefore require approximately 
12 days to complete the image combination necessary for a 
full dark energy survey. While it must be stressed that this 
is very much a "back of the envelope" estimate, and subject 
to usual additional factors due to practical realities and data 
management, the resources required are conceivably within 
reach even today (if costly). As of November 2010 this ma- 
chine was rated as the seventh fastest on earth but, thanks in 
part to the advent of inexpensive Graphical Processing Units, 
the cost and availability of such computing power must be ex- 
pected to change favorably over the next decade. Moreover, 
the survey imaging data also require considerable time simply 
to be collected by the telescope. This fact makes it unlikely 
that optimal linear image combination would represent a rate- 
determining step for a dark energy mission, even given more 
modest access to supercomputing resources. 

Finally, this discussion has assumed that the linear combi- 
nation method described in this Paper would be performed 
repeatedly and independently for each and every patch of sky. 
However, departing from this general but 'brute force' ap- 
proach there are doubtless many gains in efficiency to be made 
when processing imaging data for a dark energy survey. Not 
least, image combination efforts might be focused in regions 
where significant concentrations of light have been detected 
on a coarser grid, avoiding empty regions of the sky. Also, 
there are likely to be close similarities between the PSF and 
input pixel configurations for nearby patches: this might allow 
results from costly stages in the calculation (e.g. the eigende- 
composition of Aaij) to be reused multiple times. Further 
optimizing the implementation of the linear image combina- 
tion method would make an interesting topic for future work, 
particularly as plans for a space-based imaging survey move 
closer towards reality. 

7. CONCLUSIONS 

It has been the intention of this Paper to demonstrate a 
simple, yet general, formalism for the linear combination of 
undersampled astronomical images to produce oversampled 
output. Image reconstruction based on this formalism allows 
the explicit control of any changes to the effective image PSF 
Gi{r), and noise variance J^aa, due to the combination pro- 
cess. Such an approach is likely to be a useful addition to a 
growing set of image analysis tools being built up to face un- 
precedented challenges in astronomical inference: challenges 
arising primarily in the bid to understand dark energy. 

Despite the numerical emphasis of the approach, we have 
shown how significant computational savings may be made 
via efficient implementation, with further savings being pos- 
sible. It is planned to make the prototype Fortran 95 code 
Imcom, used to generate the results of this Paper, available 
to the public; this will be either in its current form or as part 

^ http://www.lanl.gov/roadranner/ 
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of a more developed, scaled-up software package for use with 
greater volumes of data. Even using the current prototype im- 
plementation, it has been shown that the method is a feasible 
means of processing all-sky data from a dark energy mission. 

Within the linear framework presented, we have demon- 
strated some interesting results for the WFIRST design con- 
cept adopted in Section |4] and explored in more detail in Sec- 
tion |5] First, we have found that this design allows the con- 
trolled reconstruction of fully-sampled images if five input 
dithers (in a \/5 x \/5 ideal pattern) or six input dithers (ran- 
domly offset) are available. We also found that the y/E x \/5 
ideal dither and 6 exposure random dither configurations 
proved reasonably robust in Ua despite suffering 5% bad in- 
put pixel losses in random locations. Largely, these defects 
were able to be absorbed into an increase in the noise vari- 
ance Sqq,. Whether this conclusion is robust when the bad 
pixel patterns are less random remains to be seen: the con- 
tiguous spatial extent of defects such as cosmic rays is likely 
to be a complicating factor for the control of both Ua and J^aa 
in the output, and merits further investigation. 

Another result, discussed in Section |531 was the sensitivity 
of the linear image combination method to variations in the 
pixel plate scale between exposures. Plate scale variations 
of the order ^1% were seen to significantly degrade Ua and 
Eact for 6 exposure random dither It was found that adding an 
additional Gaussian smoothing kernel of standard deviation 
(T = 0.09 arcsec (0.5 input pixels) to the desired output PSF 
r(r) mitigated this degradation significantly. The dependence 
of the strength of this effect upon input image parameters such 
as the plate scale variation, input PSF, and dither strategy will 
need to be investigated more thoroughly than was possible in 
this Paper Coping efficiently with plate scale variation clearly 
represents a non-trivial challenge when using optimal linear 
techniques to combine multiple input images. 

In the examples tested in Section |5] it was also found that 
the stationarity of the noise characterized by Y,aa (also a de- 
sirable property if Ua may be kept low) seems to be best pre- 
served for the ideal dither patterns, even while Ua may show 
strong spatial variation (e.g. Section ISTt . Conversely, setting 
limits on Ua for the random dither patterns causes significant 
spatial variation in Sqq,. This effect is also something to be in- 
vestigated further when designing dither strategies for a dark 
energy survey mission. Correlated, non-stationary noise in- 
troduces noise rectification biases that must be calculated and 
removed when performing accurate photometry and shape 
measurement. Precise calibration of such biases most likely 
proceeds best via simulations or deep (high signal-to-noise) 
training data, althou gh analytic methods exist to determine t he 
leading order terms (Kaiserll2000HBemstein & Jarvisll2002l ). 

For precisely such reasons, it is suggested that inference re- 
garding galaxy sha pes might b est occur at the level of raw 
pixel images (e.g. (Miller et all l2007t iKitching et aP l2008h . 
while image stacks or other post-processed science tools 
would be used to provide crucial ancillary information such 
as fully-sampled PSF images and precise object centroid es- 
timates. The discussion of this question is appropriate and 
timely. Whatever the eventual conclusions of such a debate, 
there is undoubted value in a method for combining multi- 



ple undersampled images to generate oversampled output that 
offers control over any changes to the PSF, while simultane- 
ously suppressing noise where possible. The formalism pre- 
sented in this Paper, being linear, is the simplest approach that 
meets these objectives. 

There remain some important unanswered questions re- 
garding the use of the technique for real data. As mentioned 
early in Section|2] the input PSF Gi{r) must be known, mod- 
eled or approximated before the use of the linear image com- 
bination formalism. Any model of Gi{r) which is sensitive 
to aliasing in images of stars from the input li risks impart- 
ing similar defects to the output image Ha- This would be 
a problem if attempting non-parametric reconstruction of the 
PSF directly from stellar images, for example. However, even 
in this case, a mitigating strategy may be found. If Gi (r) does 
not vary rapidly across the instrumental field of view then im- 
ages of different stars within a small region of a given expo- 
sure will approximate multiple images of the same PSF, but 
with a variety of centroid offsets (assuming variations due to 
star color are small, and ignoring varying flux which can be 
normalized). Fully-sampled images of an "average" star for 
each small region might thus be reconstructed by applying the 
linear image combination formalism in combination with an 
accurate estimation of the individual stellar centroids. 

In practice, any model for Gi (r) can only be an estimate of 
the true convolving kernel, and the i mpact of this uncertainty 
on estimates of key observables (e.g. lPaulin-Henriksson et al] 
2008; Rowe 2010) must be quantified in context of optimized 
linear image combination. A related question that must be 
addressed is the effect of color variation in Gi (r) for images 
observed using broad-band spectral filters. Furthermore, and 
beyond its importance when modelling G,; (r), the need for ac- 
curate stellar centroid estimation also arises when considering 
astrometric registration for each of the individual exposures 
that make up li. This registration is needed before the images 
may be combined, and any model of r, can only approximate 
the true astrometric solution. The impact of realistic levels 
of uncertainty in upon the quality of the output image Ha 
must be investigated. Work to explore some of these impor- 
tant questions is already underway in the laboratory, where 
they are equally relevant in the precise characterization of de- 
tector technology for a dark energy mission. Happily, the for- 
malism presented here provides a useful framework for ex- 
ploring these issues. 
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APPENDIX 

ASYMPTOTIC BEHAVIOR OF OUTPUT IMAGE PROPERTIES 

The solution to the system of equations derived by minimizing Wa in equation (IT6] l. which is given by equation ( l20b . may be 
written succinctly as 

T = -i(A„ + KN)-^B„. (Al) 
Using this result with equation ^ we may write the noise variance in each output pixel as 

^aa —^^TaiTaj Nij (A2) 

ij 

= iBj(A„ + KN)"'N(A„ + KN)"'B„. (A3) 
The leakage objective Ua at each pixel may similarly be written 

C/„ = ^Bl (A„ + kN)-' A„ (A„ + kN)-^ - ^B^ (A, + kN)-' B„ + C„. (A4) 
Noting the simple k dependence in the expressions above, we see immediately that 

— ^ = --B^ (A„ + kN)-^ N (A„ + kN)-' N (A„ + kN)-' B„ < 0, (A5) 

and 

^ = \bI (a„ + ^ny^ N [i„ - (A„ + uny^ aJ (a„ + n^y^ b„ > o, (A6) 

where I„ is the n x n identity matrix. 

These results can be interpreted in terms of the asymptotic behavior of the output image properties across the full range of k. 
Taking first the limit k oo,we see that 

lim{S„a}^0+, lim|^|^l^O-. (A7) 



Ok 



As K becomes large the system ten ds co nvergently to an output image of zero noise variance (or covariance). This result can in 
fact be trivially seen from equation (lAll i. as T in the limit k —?' oo. The output image Ha likewise tends to zero everywhere, 
and the leakage and leakage objective tend as 

lim {L„(r)}^-r(r), lim ^ C„, lim(^1^0+. (A8) 



Ok 



It follows hence that values of the leakage objective Ua are most naturally quoted in units of Cq, converging upon unity in the 
limit of least image fidelity. 
Turning to opposing limit k 0, we find that 

lim{Ua}^Ca-jBlAa'Ba, lim|^1^0+. (A9) 

K-i-0 4 K^o I OK ) 

While we can see, therefore, that the system converges stably to this minimum value there is no guarantee that J/q, — > as k — > 0, 
and therefore no guarantee of an unbiased linear solution for Tai- The noise variance exhibits the following hmiting behavior 

lim {E„J ^ iB^A-iNA-iB„, lim j ^ -iB^A-^NA-^NA-iB^. (AlO) 

Depending on Aaij and Bai, these variances may be extremely large resulting in noisy output images. In practice a compromise 
between noise and fidelity will be desired. 
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